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ALLMAT: A TSS/360 FORTRAN IV SUBROUTINE FOR EIGENVALUES 
AND EIGENVECTORS OF A GENERAL COMPLEX MATRIX 

by Gale Fair 
Lewis Research Center 

SUMMARY 

A subroutine is described and listed that computes the eigenvalues and eigenvectors 
of a general (non-Hermitian) complex matrix. The program, ALLMAT, uses the com- 
plex QR algorithm to compute eigenvalues and inverse iteration to compute eigenvectors. 
The user has the option of computing only the eigenvalues, if desired. An entry point 
EVDATA is available to provide the user with timing and accuracy information, as well 
as the number of iterations necessary for each eigenvalue and eigenvector. 


INTRODUCTION 

Many areas of physics, mathematics, statistics, and engineering require the eigen- 
values and eigenvectors of square matrices. This area of numerical analysis, some- 
times called the algebraic eigenvalue problem, holds a place that is as important as the 
more familiar areas such as numerical integration, curve fitting, and numerical in- 
tegration of differential equations. A general library of subroutines for a computer 
installation is commonly limited to one such program, and quite often this one sub- 
routine is of only limited applicability. 

The ideal subroutine for the algebraic eigenvalue problem should have many fea- 
tures: it should be fast and accurate; it should give a matrix of eigenvectors that are 
linearly independent; it should be capable of computing only eigenvalues at a corre- 
sponding increase in speed; it should have minimal storage requirements; and it must 
be able to treat all matrices, real or complex, symmetric or nonsymmetric, regardless 
of the condition of the matrix. Unfortunately, probably no such procedure exists. The 
usefulness of any subroutine may be judged on the basis of how many of these criteria 
are fulfilled, as balanced against the needs of the individual user. 


There are many different techniques for diagonalizing a square matrix. The 
treatise by Wilkinson (ref. 1) is evidence for this. Reference 1 describes the state-of- 
the-art for the algebraic eigenvalue problem as of 1965. Typically, one chooses a 
particular method because he believes that his matrix has some feature that requires 
special handling, because a subroutine is conveniently at hand, or because he knows only 
that one method. Two of the most commonly used procedures are the power method and 
the Jacobi transformation. 

The power method is a special-purpose procedure that computes the largest eigen- 
value of a matrix by the formation of a sequence of powers of the matrix acting upon an 
arbitrary vector. This procedure is useful for the computation of a few eigenvalues 
(the largest in magnitude) and their eigenvectors. The computation of a full set of 
eigenvalues and eigenvectors is both time consuming and inaccurate. 

The Jacobi transformation, as it applies to a complex Hermitian matrix, consists 
of a sequence of unitary transformations that diagonalize 2x2 submatrices of the full 
matrix. This procedure generates the eigenvectors along with the eigenvalues and is 
particularly useful when the eigenvectors are required to be orthogonal to a high degree 
of accuracy. The limitations of the Jacobi transformation are that the accuracy of the 
eigenvectors is usually limited, and as yet no extension to nonsymmetric or non- 
Hermitian matrices has been made. The most common computer library subroutine 
for the algebraic eigenvalue problem is a real-symmetric version of the Jacobi method 
(ref. 2). 

For a general (i.e. , nonsymmetric or non-Hermitian) matrix, two procedures have 
been derived to compute the eigenvalues and eigenvectors, respectively. The input 
matrix is reduced to a Hessenberg form (ref. 1), and the QR transformation of Francis 
(refs. 3 and 4) is used to compute the eigenvalues. With a knowledge of the eigenvalues, 
the Wielandt inverse iteration method (ref. 1) generates the eigenvectors. The QR 
transformation and inverse iteration appear to be the best currently available for their 
respective tasks (ref. 1) in terms of accuracy and speed. This combined procedure has 
been coded at the Oak Ridge National Laboratory (ref. 5) for an IBM 360/50 using the 
H-level FORTRAN compiler and COMPLEX*16 arithmetic. This program was used as 
the basis for the subroutine to be described in this report. 

In order to make the subroutine as general as possible, some modifications and 
additions were made to the ORNL program, as follows. The fact that for some matrices 
the Hessenberg form may be decomposed into disjoint submatrices is incorporated in 
both the QR transformation and the inverse iteration to reduce computational time. A 
perturbation method is used to obtain linearly independent eigenvectors when eigenvalues 
are either degenerate or very nearly the same value (meaning that the matrix itself may 
be ill-conditioned (ref. 1)). An auxiliary entry point is provided to give the user in- 
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formation about the number of iterations required, timing data (measured as central 
processor or CPU time elapsed in the computation), and error data for the resulting 
eigenvalues and eigenvectors. Finally, a flag has been provided to allow the user to 
compute only the eigenvalues, with the use of the QR transformation. The relative con- 
tribution made by the present work is seen from the observation that approximately 
60 percent of the coding of the current form of the subroutine ALLMAT is the ORNL 
coding while the remaining 40 percent is new. 

The end result of the work described here is a subroutine for the IBM/360 to com- 
pute the eigenvalues and eigenvectors of a square matrix. Certainly, this is not meant 
to be the final work in such procedures, the algebraic eigenvalue problem is an area 
of extensive research in numerical analysis. On the other hand, this subroutine does 
satisfy most of the criteria mentioned earlier for the ideal subroutine, at least to some 
degree. The criterion that is least satisfied is minimal storage. Because ALLMAT is 
written with COMPLEX*16 arithmetic and has some large scratch-pad arrays, the sub- 
routine uses a large amount of storage. On a TSS/360 system this storage requirement 
is not a basic limitation on the subroutine, but it does imply that the CPU time is af- 
fected. 

This brief mention of storage requirements is an opportunity to interpose a slight 
warning to the prospective user of ALLMAT. If only a small number of the (largest) 
eigenvalues of a matrix are desired, the power method is more efficient than ALLMAT. 
For a real, symmetric matrix, a problem that requires eigenvectors along with the 
eigenvalues would be better suited to a real Jacobi subroutine. On the other hand, for 
the computation of eigenvalues alone, or for the eigenvalues and eigenvectors of a real, 
nonsymmetric matrix or for a complex matrix, ALLMAT seems to be the best choice, 
at this time. 

The next section of this report describes schematically the construction of the sub- 
routine ALLMAT. This includes the information necessary for a programmer to use 
ALLMAT. Also included are brief descriptions of the mathematical procedures used in 
ALLMAT. The following section discusses the special features that have been incorpo- 
rated in ALLMAT, including a description of the subsidiary ENTRY EVDATA that pro- 
vides timing and accuracy information for the user. Finally, a number of test matrices 
are used as examples for ALLMAT. These examples give an indication of running 
times and accuracy obtainable with the program, even with some ill-conditioned input 
matrices. A FORTRAN listing of ALLMAT is given in the appendix. 

This report is intended to be used as a user's manual for the subroutine ALLMAT, 
and as such it described the call vector for the subroutine and the rules for usage. In 
addition, enough information is provided the prospective user to allow an intelligent 
application of this program to his particular problem. The prospective user should not 
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apply this program to his problem without some understanding of the numerical methods 
involved and of the construction of the subroutine. 


GENERAL CONSTRUCTION 
Usage 

The information to be discussed in this section is aimed at explaining the program 
as a FORTRAN subroutine, along with a description of the ENTRY EVDATA. 

The user’s access is through the statement (see the appendix for the complete 
FORTRAN listing of the subroutine): 

CALL ALLMAT (AA, LAMBDA, M, MM, EVECT, NCAL) 


where 

AA input COMPLEX*16 matrix, of dimension M ^ MM. Upon return from 

ALLMAT, i th column of AA is i th eigenvector, corresponding to i th 
eigenvalue. 

LAMBDA COMPLEX* 16 vector of length M that contains eigenvalues upon return 
from ALLMAT. 

M actual dimension of input matrix AA. 

MM dimension of AA as it appears in a dimension statement in the calling pro- 

gram. MM is the upper bound for the size of matrices used. As 
ALLMAT is currently written, MM must be no greater than 50. 

EVECT a logical switch. If EVECT = .TRUE. , the eigenvectors of AA are calcu- 
lated, and returned in the matrix AA. If EVECT = .FALSE. , no eigen- 
vectors are calculated and AA contains no useful information upon return 
from ALLMAT. 


NCAL number of eigenvalues successfully computed by ALLMAT. If NCAL < M 

some attempts of the QR transformation did not converge within 10' iter- 
ations. The value of the element of LAMBDA that corresponds to this 
eigenvalue has been set to zero by ALLMAT. 


In addition to the primary entry point, a secondary ENTRY EVDATA is available to 
give the user information on the CPU time taken for the eigenvalue and eigenvector 
procedures. Also available are the number of QR iterations required for each eigen- 


4 


I 


value, the number of inverse iterations required for each eigenvector, and the Euclidean 
norms of the residual vectors. A more complete description of these quantities is given 
later. The usage for this optional entry point is 

CALL EVDATA (ITS,KTS,NCO, MCO,RNORM) 

where 

ITS elapsed time for QR transformation for eigenvalues, including time to re- 

duce to upper Hessenberg form. ITS is an integer, in microminutes. 

KTS elapsed time for inverse iteration for eigenvectors. Does not include time 

represented by ITS. Also an integer in microminutes. 

th 

NCO an integer vector of dimension MM that has as its 1 11 element the number 

of QR iterations for the i^ 1 eigenvalue. NCO (i) s 10. If NCO (i) = 0, 
this eigenvalue was obtained along with another, no separate QR iteration 
was required. If NCO (i) < 0, no convergence was obtained for this eigen- 
value within ten QR iterations. 

th 

MCO integer vector of dimension MM that has as its i n element the number of 

th 

inverse iterations necessary to obtain the i m eigenvector. MCO (i) ^ 10. 

RNORM REAL*8 vector of the norms of the residual vectors of AA. See section 
SPECIAL FEATURES OF ALLMAT for a more complete description. 
RNORM also has a dimension MM. 

As an example of the usage of ALLMAT, consider a 6x6 complex matrix AA that 
is to be diagonalized. Let us assume that the TYPE statement in the calling program 
that specifies the dimensions of AA and LAMBDA has the form 

COMPLEX*16 AA (10, 10), LAMBDA(IO) 

The arrays have been overdimensioned for more generality. Let us further assume 
that eigenvectors are desired from ALLMAT, so that EVECT has been assigned a value 
. TRUE. . Then the call to ALL MAT is 

CALL ALLMAT (AA, LAMBDA, 6, 10, EVECT, NCAL) 

Upon return from ALLMAT the integer variable NCAL contains the number of eigen- 
values that have been successfully computed by ALLMAT. The i^ 1 column of AA (I. E. 
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xt_ 

AA (1, 1) to AA (6, 1)) contains the i 1 eigenvector, corresponding to the eigenvalue 
LAMBDA (1). 

If the timing and error information provided by EVDATA are desired by the user, 
then the statement 


CALL EVDATA (ITS,KTS,NCO,MCO,RNORM) 

is used, where NCO, MCO, and RNORM have been dimensioned at least six in the calling 
program. The conversion from ITS or KTS (in microminutes) to milliseconds is ob- 
tained by multiplying either integer by 0. 06 and assigning the result to a floating-point 
variable. 


QR TRANSFORMATION 

The basis of the QR transformation is a theorem by Francis that states any non- 
singular matrix A has a unique decomposition into the product of a unitary matrix Q 
and an upper triangular matrix R (ref. 3), or 


A = QR 


The QR algorithm consists of forming a sequence of matrices similar to A (=A^ ) 
such that 


A (K) “ Q (K) r (k) 


and then 


A (K+1) " R (K) Q (K) 


where A^ is the form of the matrix after the K 


th 


decomposition. Francis (ref. 3) 
shows that this sequence of matrices has as its limit an upper triangular matrix, the 
diagonal elements of which are the eigenvalues of the original matrix A. Furthermore, 
even if the original matrix is singular, the algorithm still gives convergence to a unique 
triangular matrix, even though some of the intermediate Q and R may not be unique. 

A full description of the QR transformation is certainly not relevant to this report. 
A detailed discussion of the convergence properties and the error analysis of the QR 
algorithm is given in references 1, 3, and 4. It is sufficient to note for our purpose 
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that the QR algorithm is an extremely stable, rapidly converging procedure to calculate 
the eigenvalues of a general matrix (ref. 1). The version of the QR transformation 
that is part of ALLMAT, one that includes origin shifts to accelerate convergence, is 
powerful enough to satisfy nearly all of the needs of the average user. 

There is one unusual feature of the standard way in which the QR algorithm is em- 
ployed that the prospective user should be aware of. A preliminary step in any im- 
plementation of the QR transformation is the reduction of the input matrix to Hessenberg 
form. An upper Hessenberg form (i. e. , A.. = 0 if i > j + 1) is used in ALLMAT. The 
reduction is accomplished by a sequence of elementary transformations (ref. 1). The 
elements of these elementary transformations are stored in the unused portion of A (the 
lower subtriangle of A) and in the integer vector JNT. This information is used at the 
end of the inverse iteration to recover the eigenvectors of the original matrix from the 
eigenvectors of the Hessenberg matrix. The point of caution for the user is that the 
working matrix for the subroutine is the Hessenberg form, which in general bears no 
simple relation to the input matrix. Thus, if the user attempts to debug this subroutine 
at an intermediate stage, the relation between the Hessenberg form and the original form 
must ge kept in mind. 

The advantage of using the Hessenberg form is apparent in the time needed to com- 
plete the computation of the eigenvalues. Most methods that operate on the entire input 
matrix, such as the Jacobi method, require a number of operations that is approxi- 

O 

mately 30N (ref. 2), where N is the order of the matrix. The reduction to Hessenberg 

3 

form is a one-pass operation and requires o;N operations, where a is of order unity. 

The QR algorithm applied to the Hessenberg form only requires something of the order of 
o 

N operations. One interesting result of this is the observation (ref. 5) that under many 
conditions the QR transform produces eigenvalues in less time than the Jacobi trans- 
formation. 


Inverse Iteration 

The basis of the inverse iteration procedure is the observation that, if X is an 

eigenvalue of the matrix A, the quantity (A - XI), where I is the unit matrix, will be 

singular. Thus, if X is a good approximation to an eigenvalue of A, the matrix 

(A - XI) -1 may be iterated to obtain an approximation Y to the eigenvector X. The 

th 

iteration process is carried out until after the K n iteration the norm of the iterated 
vector, (A - XI) - * is greater than some preselected value (see the appendix). This 
procedure is equivalent to the power method, but in inverse powers of the matrix 
(A - XI). The speed with which this iteration produces an eigenvector depends on the 
accuracy of the estimate for the eigenvalue, but rarely does this procedure, combined 
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with the QR algorithm, require more than 2 iterations to produce eigenvectors to at 
least six or seven place accuracy. Again, the interested user is referred to Wilkinson 
(ref. 1) for a complete description of the method and the error analysis. 


SPECIAL FEATURES OF ALLMAT 

As mentioned in the introduction, the basic elements of ALLMAT, the reduction to 
Hessenberg form, the QR transformation , and the inverse iteration, are taken from an 
ORNL subroutine (ref. 5). There are several features that have been added to this basic 
program to either add effectiveness to the program or provide timing and accuracy in- 
formation to the user. These special features will be discussed in this section, more 
or less in the order that they appear in the program. 


Decomposed Hessenberg Form 

The reduction of the original matrix to Hessenberg form is a procedure that de- 
creases the number of operations necessary for the QR algorithm. In a large number 
of cases the nature of the Hessenberg form allows further simplifications. To illustrate 
this, sketch (a) shows an upper Hessenberg matrix, of order N. The X’s in the sketch 
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indicate matrix elements, generally nonzero, whose values are unimportant. Now let 
one of the subdiagonal elements vanish, for example A(R, R - 1) = 0. Then the 
Hessenberg matrix may be decomposed into four submatrices as shown in sketch (b). 



(b) 

The submatrices B and D are upper Hessenberg matrices of order R - 1 and 
N -R + 1, respectively. Submatrix C is a nonzero matrix with N - R + 1 columns 
and R - 1 rows. The remaining submatrix of this partition of A is entirely filled 
with zeroes. 

The result of this decomposition is that the problem of finding the eigenvalues of 
B and D becomes entire disjoint; that is, the eigenvalues of B and D, collectively, 
are the eigenvalues of A. The submatrix C plays no part in the eigenvalue problem. 
Thus, instead of the solution of a single matrix of order N, the problem has been re- 
duced to the solution of two matrices, of order N - R + 1 and R - 1. Since 

> (N - R + l)^ + (R - 1)^ for N^3 and R that is not trivial, this decomposition 
implies a significant reduction in the total number of operations in the QR transforma- 
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tion. For many input matrices, particularly those matrices that are sparse, a number 
of such decompositions may be performed and the gain in machine time is important. 

This decomposition is not as important for the calculation of the eigenvectors, 
although some improvement is made. The eigenvectors corresponding to eigenvalues 
of D (see sketch (b)) depend upon the submatrices B and C, so that the entire matrix 
must be used in the inverse iteration procedure. The eigenvectors corresponding to 
the eigenvalues of B, on the other hand, do not require matrices C or D, so that only 
B is used in the inverse iteration. Thus, some advantage is gained from the decom- 
position for the calculation of the eigenvectors. On the whole, though, the main ad- 
vantage of the decomposition enters in the QR transformation. 


Perturbation of Close Eigenvalues 

One difficulty with the inverse iteration method arises when two or more eigen- 
values are very nearly the same. Since every calculated eigenvalue differs from the 
’’true" eigenvalue by an amount that depends on many factors, these eigenvalues may 
not produce linearly independent eigenvectors. The way chosen to resolve this acci- 
dental degeneracy was to perturb each successive close eigenvalue by an amount small 
enough to not disturb the convergence of the iterative procedure, but large enough to 
resolve the eigenvectors into linearly independent vectors (ref. 1). The choice of the 
perturbation, EPSIL, is arbitrary and a better choice could be made for particular types 
of matrices. 

Since the existence of close but distinct eigenvalues implies that the matrix may be 
ill-conditioned (ref. 1), the accuracy of the calculated eigenvectors will be in doubt. In 
this sense, the use of a perturbation to separate the eigenvalues is an attempt to re- 
cover some useful information from a badly posed problem. Thus, for most matrices 
encountered, the existence of close but distinct eigenvalues should be rare. The occur- 
rence of multiple eigenvalues is more common. 


Multiple Eigenvalues 

The existence of a set of multiple eigenvalue is a not uncommon occurrence in phys- 
ical problems. The existence of such a set implies that there is a subspace of eigen- 
vectors that one desires the basis vectors of. In this situation the perturbation is of 
some help. If, by the process of perturbing the degenerate eigenvalues within the in- 
verse iteration process, one can obtain a set of distinct eigenvectors, even if they are 
not linearly independent, then there is a standard solution to the problem of determining 
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the 'basis vectors. For this purpose ALLMAT takes the set of distinct eigenvectors 
produced by the perturbation technique just discussed and uses a Gram-Schmidt (ref. 1) 
orthogonalization procedure to give a set of linearly independent eigenvectors. Since 
the Gram-Schmidt process involves taking the differences of nearly equal numbers in 
many cases, the accuracy of such a procedure is less than the accuracy of an inverse 
iteration vector for a distinct eigenvalue. Again, however, this represents an attempt 
to salvage as much information as one can from an undesirable situation. In practice, 
as shall be seen in the section TESTS, the results of this perturbation and orthogo- 
nalization procedure are good. 


ENTRY EVDATA 

The remaining special feature of ALLMAT is represented by the secondary entry 
point, EVDATA, as discussed in general construction. A typical user of an installation- 
supplied mathematical subroutine is usually blissfully unaware of any error consider- 
ations for his problem. Since the accuracy of any matrix eigenvalue evaluation strongly 
depends upon the properties of the input matrix, ignoring error information is equi- 
valent to shutting one’s eyes to avoid an oncoming truck. Additionally, since some 
eigenvalues and eigenvectors may in fact be absent due to nonconvergence either in QR 
or inverse iteration, the information provided by EVDATA is-important to a user. The 
use of the TSS/FORTRAN multiple-data set capability means that this information is 
readily available to the user, without so much as the disturbance of an artistic output 
format. 

The information available in EVDATA includes the number of iterations, the CPU 
time elapsed for the eigenvalue and the eigenvector computations, and an error estimate 
for each eigenvalue- eigenvector pair. The timing and counting variable provided in 
EVDATA were discussed sufficiently under usage, but the error information requires 
some further comment. 

If X and X are an exact eigenvalue and an exact eigenvector of the matrix A, then 
the vector AX - XX will be identically zero. Since neither X nor X can ever be 
computed exactly, this vector (AX - XX), called the residual vector, will be nonzero. 
The magnitude of this vector is then a measure of the error in X and X. The length of 

a vector, as used in ALLMAT, is the Euclidean norm, | |x| | = ^C|x(i)| . The vec- 

tor RNORM of EVDATA contains the norm of the residual vector for each eigenvalue- 
eigenvector pair, scaled to the Euclidean norm of the input matrix. 

The data entry point EVDATA may be used even if no eigenvectors are computed 
(i. e. , if EVECT = . FALSE. ) In this case only ITS and NCO contain meaningful values. 


11 



TESTS 


Seven matrices were chosen as examples for ALLMAT. The dimensions of these 
matrices vary from four to 19. All but two matrices are real but not symmetric, one 
of the remaining matrices is Hermitian, and the final example matrix is complex, but 
not Hermitian. Some of these matrices were chosen to illustrate ill- conditioning of one 
type or another. Since the numerical values of the eigenvectors are not of general use, 
they are not displayed. 


Matrix 1 



This real, but unsymmetric, matrix of order four has the exact eigenvalues 0. 9, 0.6, 
-0.3, and 0. In addition, the computed eigenvalue corresponding to 0. is 0. 14E-16. 
The information available from EVDATA on this test includes: 


Eigenvalue 

Number of QR 
iterations 

Number of inverse 
iterations 

RNORM 

-0.3 

8 

3 

0.65E-16 

. 14E-16 

1 

3 

. 43E-16 

.6 

1 

3 

. 16E-16 

.9 

1 

3 

. 91E-16 


The total time for the QR transformation, including the initial reduction to Hessenberg 
form was 0.053 second, and the time for the inverse iteration was 0.046 second. 
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Matrix 2 


A 2 - 


/o. 25000025 

-1.0 

-0.49999975 

-0.99999975 

-0.50000050 

0.5 

0.24999950 

0.25000025 

1.00000025 

1.25 

1.00000025 

1.25000025 

-0.4999975 

-0.25 

0.25000025 

0.50000025 


This matrix has eigenvectors identical to those of matrix 1, but has a different set of 
eigenvalues. The exact eigenvalues of Ag are given in the following table: 


Eigenvalue 

Number of QR 

Number of inverse 

RNORM 


iterations 

iterations 


1. 5 

2 

3 

0.68E-16 

. 75000075 

4 

3 

. 10E-15 

. 75 

0 

3 

. 46E-14 

75 

1 

3 

. 64E-15 


The closeness of the second and third eigenvalues hint at some error problems with the 
eigenvectors. The time for QR transformation was 0.020 second, and the time for the 
inverse iteration was 0.090 second. The difficulties anticipated from the closeness of 
the eigenvalues are evidenced in the degradation of the third value of RNORM in the 
table. 


Matrix 3 


A 3 “ 


009 

5.00101 

-8.999 

5.999 

001 

5.01101 

-8.999 

3.999 

001 

4.91101 

-8.899 

3.999 

001 

4.96101 

-8.999 

4.049 
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This matrix A g is an example of an ill conditioned matrix, in contrast with the pre- 
vious example. Here, A 2 had two nearly alike eigenvalues even though the matrix is 
not mathematically ill conditioned (ref. 1). 


Eigenvalue 

Number of QR 
iterations 

Number of inverse 
iterations 

RNORM 

0.01 

3 

3 

0.25E-16 

.01001 

3 

3 

.14E-16 

.1 

3 

3 

. 18E-16 

.05 

1 

3 

. 23E-16 


The time for QR transformation was 0.038 second; that for inverse iteration was 0.025 
second. Apparently, the ill conditioning did not effect the inverse iterations, as all 
values of RNORM are satisfactory. 


Matrix 4 



Unlike the first three test matrices, A 4 has a pair of complex eigenvalues. 


Eigenvalue 

Number of QR 
iterations 

Number of inverse 
iterations 

RNORM 

3.0929 

6 

3 

0. 78E-16 

. 1772+. 95E-16i 

6 

3 

. 24E-15 

.42295+4.39541 

5 

3 

. 1 IE- 15 

. 42295-4. 3954i 

4 

3 

. 17E-15 

15.247+. 11E-141 

1 

3 

. 56E-15 

-7. 3630 

1 

3 

. 47E-15 
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The time for QR transformation was 0. 142 second; that for inverse iteration was 0. 314 
second. The imaginary part of the sum of the eigenvalues (which should be 0. ) is 
0. 355E-14. 


Matrix 5 

This test matrix is a 19 by 19 real, unsymmetric matrix given by Francis (ref. 4) 
to demonstrate the QR transformation. The matrix is too complicated to list here, but 
the error information is informative. The time to produce the eigenvalues was 
2 seconds, and the time to calculate the 19 eigenvectors was 22 seconds. Although this 
time is large when compared with the previous examples, it is quite reasonable when 
compared with other methods (ref. 5). Even with a matrix of this order, the residual 
vectors all had norms less than l.E-16. 


Matrix 6 


0 . 

. 01375103 i 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

913751031 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

-.388213061 

0 . 

0 . 

— . 14433376 i 

0 . 

0 . 

-. 946 7409 01 ~ 

0 . 

0 . . 

38821306 i 0. 

0 . 

-.144333761 

0 . 

0 . 

-.94674090i 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

.144333761 

0 . 

0 . 

-1. 02140683 i 

0 . 

0 . 

. 03966810 r * 

0 . 

0 . . 

14433376 i 0. 

0 . 

1. 02140683 i 

0 . 

0 . 

039668101 

0 . 

0 . 

0 . 

0 . 

. 946 74090 i 

0 . 

0 . 

. 03966810 i 

0 . 

0 . 

.74705540i “ 

0 . 

0 . . 

94674090i 0. 

0 . 

-.039668101 

0 . 

0 . 

-. 74705540i 

. 0 . 


This matrix has several features that make it useful as an example. Each nonzero 
element of Ag is a purely imaginary number and, in addition, Ag is Hermitian. Thus, 
the eigenvalues of Ag are real and, since the trace of Ag vanishes, the eigenvalues 
occur in positive-negative pairs. There is a pair of degenerate eigenvalues with the 
value 0, so that the orthogonalization procedure must be used to obtain the eigenvectors. 
Finally, Ag is sufficiently sparse that the decomposition of the Hessenberg form is 
effective in reducing the time required for the computations. 
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Eigenvalue 

Number of QR 
iterations 

Number of inverse 
iterations 

RNORM 

-0.91376103 

3 

3 

0.26E-16 

.91376103 

3 

3 

. 16E-16 

-1.03872417 

7 

3 

. 31E-15 

-.38452612 

6 

3 

. 96E- 16 

.38452612 

5 

3 

. 55E-15 

1.03872417 

4 

3 

. 74E-15 

-1.53711192 

1 

3 

. 15E-14 

1.53711192 

1 

3 

. 21E-14 

0 

3 

3 

0 

0 

1 

3_ 

0 


The time for QR transformation was 1.07 second; that for inverse iteration was 
1. 50 second. The eigenvalues appear in this table in the order in which they are calcu- 
lated by the QR algorithm. Since the reduction to Hessenberg form and the use of the 
decomposed Hessenberg form rearrange the matrix, the eigenvalues are not computed 
in pairs, necessarily. This same effect caused the QR routine to take three iterations 
to compute a zero eigenvalue. The degenerate eigenvalues caused no loss of accuracy 
in the computation of the eigenvectors. Furthermore, the sum of the eigenvalues is 
purely imaginary, and has the magnitude 0. 4E-14, reflecting the zero trace of Ag. 


Matrix 7 

The final example matrix was generated from matrix by taking each element of 
this 6 by 6 real, nonsymmetric matrix and multiplying by the imaginary unit i. The 
result, Arj, is a complex non-Hermitian matrix whose eigenvalues are the eigenvalues 


Eigenvalue 

Number of QR 
iterations 

Number of inverse 
iterations 

RNORM 

3. 0929i 

6 

3 

0. 77E-16 

. 95E-16+. 1772i 

6 

3 

. 28E-15 

4. 3954+. 42295i 

5 

3 

. 10E-15 

-4. 3954+. 42295i 

4 

3 

. 16E-15 

. 44E-15-7. 3630i 

1 

3 

. 48E-15 

15. 247i 

1 

3 

. 52E-15 
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of multiplied by i. The Qr transformation time was 0. 140 second, and the inverse 
iteration time was 0. 319 second. 

A comparison of the results indicated in the preceding table with the results for 
example A^ shows that ALLMAT handles the non-Hermitian form with comparable 
speed, at no loss of accuracy in the eigenvalues and eigenvectors. The sum of the 
eigenvalues is 0. 31E-14+12. i. 

These seven examples were chosen to be representative of the application of 
ALLMAT. Some matrices (ApA^, Ag, and A^) pose no particular problems, while the 
remaining (Ag^g, and Ag) were included to demonstrate one of more special charac- 
teristics of the program. It is seen from the results given above that ALLMAT had no 
difficulty with any of these test matrices. The norm of the residual vectors is typically 
less than l.E-15, and all computed eigenvalues that were also known exactly were in 
agreement to at least 14 places. At no point did either the QR algorithm or the inverse 
iteration fail to give convergence within the allotted limit of 10 iterations. In fact, only 
once did the inverse iteration procedure require more than three iterations to satisfy 
the convergence criterion. 


CONCLUDING REMARKS 

This report is intended to be a user's guide for the prospective user of ALLMAT. 
The information presented here about the construction of ALL MAT should be considered 
a minimum for the use of this matrix eigenvector program. No program of the com- 
plexity of ALL MAT should be used without some understanding of the basic algorithms 
involved. Certainly, though, most users will apply ALLMAT without consideration of 
even the simplified discussion presented here. For these users the entry point EVDATA 
should be required usage as an indicator when ALLMAT does fail on a matrix. 

Experience has shown that two inverse iterations are usually enough to give an 
eigenvector correct to sufficient accuracy. The current version of ALLMAT, however, 
iterates until the norm of the iterated vector, (A - XI) - * X, is greater than 1.E40. If 
computing time is at a premium, this criterion can be easily changed to a test on the 
number of iterations. The current limitation on ALLMAT is to matrices of dimension 
no larger than 50. This restriction may also be changed easily. 

ALLMAT was designed to be a general purpose matrix eigenvalue and eigenvector 
subroutine. Almost any matrix, including the most general case of a complex, non- 
Hermitian matrix, is amenable to diagonalization by ALLMAT. Furthermore, timing 
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test (ref. 5) indicate that the QR transform may be preferred to the Jacobi method for 
the eigenvalues of real and symmetric matrices. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, September 16, 1970, 

129-02. 
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APPENDIX - FORTRAN LISTING 


0000100 
0000200 
0000300 
0 0 004 0 0C 
0 0 0050 0C 
0000600C 
0 00 070 0C 
000 08 00C 
0000900 
0001000 
0001100 
0001200 
0001300 
0001400 
0001500 
0001600 
0001700 100 
0001800 
0001900 
0002000 
0002100 
0002200C 
0002300C 
0 002400C 
0002500C 
0002600C 
0002700 
0002800C 
0002900C 
0003000C 
0003100C 
0003200C 
0003300C 
0003400 
0003500 
0003600 
0003700 
0003800 
0003900 
0004000 
0004100 
0004200 
0004300 


SUBROUTINE ALLMAT ( A A* LAMBDA, M,MM, EVECT, NCAL) 

-IMPLICIT REAL*8 (A-H,0-Z) 

- COMPLEX* 16 AA(MM,MM) 

(p 

IF THE USER REQU I RES^^Jl MENS I ON LARGER THAN 50, THE 
DIMENSIONS IN THE 2 LINES FOLLOWING THIS COMMENT MUST 
BE CHANGED FROM 50 TO A SIZE TO SUIT THE USER. 

■ COMPLEX* 16 A( 50, 5 0,),H(50, 50),HL(50, 5 0 ), LAMBDA (MM) 

' COMPLEX* 16 VECT(50),MULT(50),SHIFT(3),TEMP,SINT,COST,TEMP1,TEMP2 
COMPLEX* 16 E I n / CON J I 
LOGICAL EVECT, I NTH (50) 

DIMENSION N COUNT ( 50 ) , MCOIJNT ( 50 ) 

INTEGER JNT(50),R,RP1,RP2 
DO 1000 J = 1,M 
DO 1000 I = 1, M 
0 A( I , J) = AA ( I , J ) 

CALL CPUTIMC I TIM) 

ITSUM = 0 
KTSUM = 0 
CONJI = ( 0 , , -1 . ) 

THE CONSTANT EPSIL DETERMINES THE CONVERGENCE OF THE 
QR ALGORITHM, AND ALSO IS THE PERTURB AT I ON PARAMETER 
FOR THE INVERSE ITERATION. 

EPSIL = 1. OD-12 

THE CONSTANT EPSMAX DETERMINES THE CONVERGENCE OF THE 
INVERSE ITERATION. THIS NUMBER IS THE LOGARITHM OF 
NORM OF THE ITERATED EIGENVECTOR THAT IS SUFFICIENT 
FOR CONVERGENCE. 

EPSMAX = 40. 

NSTOP = M 
N = NSTOP 
NSTART = 1 
MN1 = 1 
NCAL = 0 

IF (N.NE.l) GOTO 1 
LAMBDA ( 1 ) = AC 1, 1 ) 

A( 1, 1 ) = 1.0 
GO TO 92 
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000-4400 

1 

1 COUNT = 1 

0004500 


SHIFT(l) = 0. 

0004600 


IF (N.NE.2) GOTO 4 

0004700 

2 

NS P 1 = NSTART + 1 

0004800 


TEMP = (A(NSTAPT/ NSTART )+A( NS Pl / NSP1)+CDSQRT - 

0004900 


1( (A(NSTART / NSTART) +A( NS PI, NSP1) )** 2-4. *(A(NSP1 / 

0005000 


1NSP1)*A(NSTART, NSTART ) -A ( NS PI, NSTART ) *A( NSTART, 

0005100 


INS PI ) ) ) )/2. 

0005200 


RELTEM = TEMP 

0005300 


AMGTEM = CONJ 1 *TEMP 

0005400 


IF ( RELTEM. NE.O.. OR. AMGTEM. NE.O.) GOTO 3 

0005500 


LAMBDA(NSTOP) = SHIFT(l) 

0005600 


LAMBDA (MN1) = A( NSTART, NSTART ) + A ( MSP 1, NSP1 ) +SH 1 FT ( 1 ) 

0005700 


NCOUNT ( NSTOP) = 1 COUNT 

0005800 


NCOUNT(MNl) = 1 COUNT 

0005900 


GO TO 37 

0006000 

3 

LAMBDA ( NSTOP ) = TEMP + SHIFT(l) 

0006100 


LAMBDA(MN1) = (A(NSTART, NSTART) * A ( NSP1, NS PI ) - 

0006200 


1ACNSP1, NSTART) *A( NSTART, NS PI ))/( LAMBDA ( NSTOP ) 

0006300 


2-SHI FTC 1) ) +SH 1 FTC 1) 

0006400 


NCOUNT (NSTOP) = 1 COUNT 

0006500 


NCOUNT (MN1) = 1 COUNT 

0006600 


1 COUNT = 1 

0006700 


GO TO 37 

0006800C 
0006900 C 


REDUCE MATRIX A TO HESSENBERG FORM. 

0007000C 

0007100 

4 

NM2 = N-2 

0007200 


DO 15 R = 1 , NM 2 

0007300 


R PI = R+l 

0007400 


RP2 = R+2 

0007500 


ABIG = 0. 

0007600 


JNT(R) = RP1 

0007700 


DO 5 I =RP1, N 

0007800 


RELAIR = AC 1 , R ) 

0007900 


AMGAIR = CON J 1 *A ( 1 , R ) 

0008000 


ABSSQ = RELA 1 R**2 + AMGA 1 R**2 

0008100 


IF CABSSO. LE. AB 1 G) GOTO 5 

0008200 


JNTCR) = 1 

0008300 


ABIG = ABSSQ 

0008400 

5 

CONT 1 NUE 

0008500 


INTER = JNTCR) 

0008600 


IF (ABIG.EQ.O.) GOTO 15 

0008700 


IF ( 1 NTER . EQ. RP1 ) GOTO 8 

0008800 


DO 6 1 = R , N 

0008900 


TEMP = ACRP1, I ) 

0009000 


ACRP1, 1 ) = AC 1 NTER, 1 ) 

0009100 

6 

AC 1 NTER, 1 ) = TEMP 

0009200 


DO 7 1=1, N 

0009300 


TEMP = AC 1 , RP1) 

0009400 


AC 1 ,RP1) = AC 1 , 1 NTER) 

0009500 

7 

AC 1 , 1 NTER) = TEMP 

0009600 

8 

DO 9 1 = R P 2 , N 

0009700 


MULT ( 1 ) = AC I,R)/A(RP1,R) 

0009800 

9 

AC 1 ,R) = MULT ( 1 ) 

0009900 


DO 11 1=1, RP1 
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0010000 


TEMP = 0. 


0010100 


DO 10 J=RP2,N 


0010200 

10 

TEMP = TEMP + A( 1 , J)*MULTCJ) 


0010300 

11 

A( 1 , RP1) = A ( 1 ,RP1) + TEMP 


0010400 


DO 13 1 =RP2, N 


0010500 


TEMP = 0. 


0010600 


DO 12 J=RP2, N 


0010700 

12 

TEMP = TEMP + AC 1 , J ) *MIJLT ( J ) 


0010800 

13 

ACI,RP1) = A( 1 , RP1 ) +TEMP-MU LT ( 1 )*A(RP1 / RP1) 


0010900 


DO 14 1 =RP2/ N 


0011000 


DO 14 J=RP2, N 


0011100 

14 

AC 1 , J ) = AC 1 , J ) - MULT C 1 ) *ACRP1/ J) 


0011200 

15 

CONTINUE 


0011300C 




00 1 14 00C 


CALCULATE EPSILON. 


0011500C 




0011600 


EPS = 0. 


0011700 


DO 16 1 = 1, N 


0011800 

16 

EPS = EPS +CDABS C A C 1, 1 ) ) 


0011900 


DO 18 1=2, N 


0012000 


SUM = 0. 


0012100 


1 Ml = 1 - 1 


0012200 


DO 17 J = 1 M 1 , N 


0012300 

17 

SUM = SUM + CDABS C A C 1 , J ) ) 


0012400 

18 

1 F C SUM . GT . EPS ) EPS = SUM 


0012500 


EPS = 'BSQRTCRf LOATCN) )*EPS*1. D-20 


0012600 


IF CEPS.EQ.O.) EPS=l.D-20 


0012700 


EPSI L = DMAX1CEPS, EPS 1 L) 


0012800C 


V' 


0012900C 


SAVE THE HESSENBFRG FORM IN THE ARRAY H. 


0013000C 




0013100 

20 

DO 10 1 = 1, N 


0013200 


DO 19 J=1,N 


0013300 

19 

H( 1 , J) = AC 1 , J) 


0013400 


NSM 1 = NSTOP - 1 


0013500 


IF CNSM1.NE.0) GOTO 100 


0013600 


R = 1 


0013700C 




0013800C 


START SCANNING FOR ZEROES IN THE SUB-DIAGONAL. 

THIS 

0 0 1 39 0 OC 


DEFINES THE SUB-BLOCKS OF THE DECOMPOSED HESSENBFRG 

0 01 40 OOC 


FORM. 


0 0 1 4 1 0 OC 




0014200 


GO TO 102 


0014300 

100 

DO 101 1=1, NSM1 


0014400 


R = NSTOP -1+1 


0014500 


RELAM1 = AC R, R-l) 


0014600 


AMGAM1 = CONJ 1 *ACR, R-l) 


0014700 


IF C C-QABS C RELAM1) +1DABS C AMGAM1 ) ) . LE . EPS 1 L ) GOTO 

102 

0014800 

101 

CONTINUE 


0014900 


R = 1 


0015000 

102 

NSTART = R 


0015100C 




0015200C 


NSTART AND NSTOP ARE THE INDICES OF THE BEGINNING AND 

0015300C 


END OF A DECOMPOSED HESSENBERG BLOCK. 


0015400C 




0015500 


NS = NSTOP - NSTART + 1 



0015600 
0015700 
0015800 
0015900 
0016000 
0016100 
0016200 
0016300 
0016400 
0016500 
0016600 
0016700 
0016800 
0016900 
0017000 
0017100 
0017200 
0017300 
0017400 
0017500 
0017600 
0017700 
0017800 
0017900C 
0018000C 
0018100C 
0018200 
0018300 
0018400 
0018500 
0018600 
0018700 
0018800 
0018900 
0019000 
0019100 
0019200 
0 0 1 °3 00 
0019400 
0019500 
0019600 
0019700 
0019800 
0019900 
0020000 
0020100 
0020200 
0020300 
0020400 
0020500 
0020600 
0020700 
0020800 
0020900C 
0021000C 
0021100C 


NC = NS 

MN1 = NSTOP + NSTART - N 
103 IF (NS.NE.l) GOTO 21 

LAMBDA(MNl) = A( NSTART, NSTART ) + SHIFT(l) 

NCOUNT(MNl) = I COUNT 
GO TO 37 

21 I F ( NS . EQ. 2 ) GOTO 2 
/ 22 RELANN = A(N,N) 

V AMGANN = CONJ I *A ( N, N) 

-RLNNM1 = A( N, N- 1 ) 

-AMNNM1 = CONJ I *A ( N, N-l ) 

RLNON1 = ACM, N-l) /A (N, N) 

AMNDN1 = CONJ I *(A(N,N-1)/A(N,N)) 

IF (RELANN. NE . 0 .. OR .AMGANN . NE . 0 . ) - 
* 1 1 F (0rABS(RLNONl) + .DABS ( AMNDN 1 ) - 1 . D- 18 ) 24, 24, 23 

23 IF ( DABS ( RLNNM1 )+HABS ( AMNNM1 ) . GE . EPS ) GOTO 25 

24 LAMBDA(MNl) = ACM, N) + SHI FTC 1) 

NCOUNT (MN1) = I COUNT 

I COUNT = 1 
N = M - 1 
NS = NS - 1 , 

MN1 = MN1 + 1 
GO TO 21 

DETERMINE SHIFT 

25 SH I FTC 2 ) = ( AC N-l, N-l ) +A ( N, N) +CDSQRT ( C A CN- 1, - 
1N-1)+A(N,N) )**2-4.*(A(N,N)*A(M-l,N-l)-A(N, N-l) - 
2*A(N-l,N))))/2. 

RELSHF = SHIFTC2) 

AMGSHF = CONJ I *SH I FTC 2) 

IF (RELSHF. NE.O.. OR. AMGSHF. NE.O.) GOTO 26 
SH I FTC 3 ) = A(N-1,N-1)+A(N,N) 

GO TO 27 

26 SHI FT(3)=CA(N,N)*A(M-1,N-1)-A(N,N-1)*ACN-1,N) )/SHIFTC2) 

27 IF(CfcABS(SHIFT(2)-A(N,N)).LT.CbABS(SH|FT(3) - 
1-A(N,N) )) GO TO 28 

INDEX = 3 
GO TO 29 

28 INDEX = 2 

29 I F(CDABS(A(N-1,N~2 ) ) .GE. EPS) GOTO 30 
LAMBDA (MN1) = SH I FT ( 2 ) + SHIFT(l) 

LAMBDACMNl+l) = SHIFTC3) + SHIFT(l) 

NCOUNT (MN1) = I COUNT 

NCOUNT (MN1+1) = 0 
I COUNT = 1 
N = N - 2 
NS = NS - 2 
MN1 = MN1 + 2 
GO TO 103 

30 SH I FTC 1 ) = SHIFT(l) + SHIFT(INQEX) 

DO 31 I =NSTART, N 

31 AC I , I ) = AC I , I ) - SH I FT ( I NDEX ) 

PERFORM GIVENS ROTATIONS, QR ITERATES. 
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0021200 

0021300 

0021400 

0021500 

0021600 

0021700 

0021800 

0021900 

0022000 

0022100 

0022200 

0022300 

0022400 

0022500 

0022600 

0022700 

0022800 

0022900 

0023000 

0023100 

0023200 

0023300 

0023400 

0023500 

0023600 

0023700 

0023800 

0023900 

0024000 

0024100 

0024200 

0024300 

0024400 

0024500 

0024600 

0024700 

0024800C 

0024900C 

0025000C 

0025100 

0025200 

0025300 

0025400 

0025500 

0025600 

0025700 

0025800 

0025900 

0026000 

0026100 

0026200 

0026300 

0026400 

0026500 

0026600 

0026700 


IF ( I COUNT. LE. 10) GOTO 32 
NCOUNT(MNl) = -ICOUNT 
NC * NC - NS 
GO TO 37 

32 NM1 * N - 1 

TEMPI = A ( NS TART, NS TART ) 

TEMP2 = A(NSTART+1,NSTART) 

DO 36 RESTART, NM1 
NN = R 
RP1 = R + 1 
RELTM1 = TEMPI 
AMGTM1 = CONJ I *TEMP1 
RELTM2 * TEMP2 
AMGTM2 = CONJ I *TEMP2 

RHO «'«PSQRT( RELTM1**2+AMGTM1'**2+RELTM2**2 + AMGTM2**2 ) 

IF (RHO.EQ.O.) GOTO 36 
COST » TEMP1/RHO. 

SINT » TEMP2/RHO 

INDEX « MAX0(NN-1,NSTART) 

DO 33 I = I NDEX, N 

TEMP = l)CONJG(COST)*A(NN, I )+lNCONJG(SINT)*A(RPl, I ) 

A( RP1, I ) = -S I NT *A (NN, I ) + COST *A(RP1, I ) 

33 A ( NN, I ) = TEMP 
TEMPI = A(RP1,RP1) 

TEMP 2 = A ( NN + 2, R+ 1 ) 

DO 34 I =NSTART, R 

TEMP= COST * A ( I, NN)+S I NT*A ( I ,RP1 ) 

A( I , RP1 ) = -DCONJGC S I NT) *A( I , NN ) +DCON JG (COST ) *A ( I , R PI ) 

34 A( I ,NN> = TEMP 

INDEX = M I NO ( NN+ 2, N ) 

DO 35 I =RP1, I NDEX 
A ( I , NN) = S I NT*A ( ! , RP1 ) 

35 A( I , RP1 ) -DCONJG (COST) *A( I , RP1 ) 

36 CONTINUE 

! COUNT = I COUNT + 1 
GO TO 22 

CALCULATE VECTORS. 

37 IF ( . NOT . EVECT ) GOTO 64 
CALL CPUTIM(JTIM) 

ITSUM = ITSUM + (JTIM - I T I M ) 

IF (NC.EQ.O) GOTO 64 
NPNCAL = NS TART + NC - 1 
N » NSTOP 

NS = NSTOP - NSTART + 1 
NM1 = N - 1 
IF (N.NE.2) GO TO 38 

EPS = &MAXK CRABS ( LAMBDAC 1) ) , CE^ABS ( LAMBDA ( 2 ) ) ) *1. D-16 
IF (EPS.EQ.O.) EPS = E PS I L 
H ( 1, 1 ) = A ( 1 , 1 ) 

H ( 2, 1 ) = A ( 2, 1 ) 

H ( 1, 2 ) - A ( 1, 2 ) 

H ( 2, 2 ) = A C 2, 2 ) 

38 DO 63 L=NSTART, NPNCAL 
ABIG = 0. 
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0026800 


EIG = LAMBDA(L) 

0026900 


IF (L.EQ.NSTART) GOTO 40 

0027000 


LM1 = L - 1 

0027100 


RELEIG = EIG 

0027200 


AMGEIG = CONJ 1 *E 1 G 

0027300 


DO 39 1 =NSTART, LM1 

0027400 


RELAM 1 = LAMBDA ( 1 ) 

0027500 


AMGAMI = CONJ 1 *LAMBDA( 1 ) 

0027600 


IF (DABSCRELEI G- RE LAM 1 ) . GT. EPS 1 L ) GOTO 39 

0027700 


IF ( DABS (AMGE 1 G-AMGAM 1 ) .GT. EPS 1 L) GOTO 39 

0027800 


EIG = EIG + CONJ I * E PS 1 L 

0027900 

39 

CONTINUE 

0028000 

40 

DO 42 1 =1, N 

0028100 


DO 41 J = 1,N 

0028200 

41 

HL ( J, 1 ) = H(J, 1 ) 

0028300 

42 

HL( 1, 1 ) = HL( 1 / 1 ) - EIG 

0028400 


DO 46 1 =1/ NM1 

0028500 


MULT ( 1 ) = 0. 

0028600 


1 NTH ( 1 ) = .FALSE. 

0028700 


1 PI = 1 + 1 

0028800 


IF (CDABS(HL(I+1, l)).LE.CDABS(HL(l, 1))) GO TO 44 

0028900 


INTH(I) = .TRUE. 

0029000 


DO 43 J= 1 , N 

0029100 


TEMP = HL( 1 +1, J) 

0029200 


HL(I+1,J) = HL( 1 , J) 

0029300 

43 

HL ( 1 / J) = TEMP 

0029400 

44 

RELH! ! = HL ( 1 / 1 ) 

0029500 


AMGHI 1 = CONJ l*HL( 1,1) 

0029600 


IF (RELHI 1 .EQ.O.. AND. AMGHI 1 .EQ.O. ) GOTO 46 

0029700 


MULT ( 1 ) = -HL( 1+1, 1 )/HL(l, 1 ) 

0029800 


DO 45 J=IP1,N 

0029900 

45 

HLU+l, J)=HL( 1 + 1, J) + MULT ( 1 ) *HL ( 1 , J ) 

0030000 

46 

CONTINUE 

0030100 


DO 48 1=1, N 

0030200 

48 

VECT ( 1 ) = 1. 

0030300 


IF (NSTOP.EQ.M) GOTO 110 

0030400 


NSTP1 = NSTOP + 1 

0030500 


DO 47 1 =NSTP1,M 

0030600 

47 

VECT ( 1 ) = 0. 

0030700 

110 

1 COUNT = 1 

0030800 

49 

RELHNN = HL(N, N) 

0030900 


AMGHNN = CONJ 1 *HL (N, N) 

0031000 


IF (RELHNN. EQ.O. .AND. AMGHNN. EQ.O.) HL(N,N)=EPS 

0031100 


VECT(N) = VECT (N)/HL(N,N) 

0031200 


DO 51 1=1, NM1 

0031300 


K = N— t 

0031400 


DO 50 J=K, NM1 

0031500 

50 

VECT(K) = VECT(K) - HL ( K, J+ 1 ) *VECT( J+ 1 ) 

0031600 


RELHKK = HL ( K, K ) 

0031700 


AMGHKK = CONJ 1 *HL(K,K) 

0031800 


IF (RELHKK. EQ.O. .AND. AMGHKK. EQ. 0. ) HL(K,K)=EPS 

0031900 

51 

VECT(K) = VECT (K) /HL (K, K) 

0032000 


BIG = 0. 

0032100 


DO 52 1=1, N 

0032200 


RELVEC = VECT ( 1 ) 

0032300 


AMGVEC = CONJ 1 *VECT ( 1 ) 
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0032400 SUM = DABS (RELVEC)+DABS(AMGVEC) 

0032500 IF (SUM. LE. BIG) GOTO 52 

0032600 BIG = SUM 

0032700 11=1 

0032800 RELV = RELVEC 

0032900 AMGV = AMGVEC 

0033000 52 CONTINUE 

0033100 IF (BIG.EQ.O.) GOTO 155 

0033200 IF (AMGV.EQ.O.) GOTO 135 

0033300 IF ( DABS (AMGV ).GT. DABS (RELV)) GOTO 125 

0033400 RAT = AMGV/RELV 

0033500 DEN = RELV + RAT *AMGV 

0033600 DO 120 1=1, N 

0033700 IF (I .EQ. I I ) GOTO 120 

0033800 RELVEC = VECT(I) 

0033900 * AMGVEC = CON J I *VECT ( I ) 

0034000 RELVC = (RELVEC + RAT* AMGV EC )/DEN 

0034100 AMGVC = (AMGVEC - RAT*RELVEC) /DEM 

0034200 VECT(I) = DCMPLX (RELVC, AMGVC) 

0034300 120 CONTINUE 

0034400 VECT ( I I ) = 1. 

0034500 GO TO 150 

0034600 125 RAT = RELV/AMGV 

0034700 DEN = AMGV + RAT*RELV 

0034800 DO 130 1=1, N 

0034900 IF ( I . EQ. I I ) GOTO 130 

0035000 RELVEC = VECT ( I ) 

0035100 AMGVEC = CON J I *VECT ( I ) 

0035200 RELVC = (AMGVEC + RAT*RELVEC ) /DEN 

0035300 AMGVC = ( RAT *AMGVEC - RELVEO/DEN 

0035400 VECT ( I ) = DCMPLX ( RELVC, AMGVC ) 

0035500 130 CONTINUE 

0035600 VECT ( II) = 1. 

0035700 GO TO 150 

0035800 135 DO 53 1=1, N 

0035900 53 VECT(I) = VECTCD/BIG 

0036000 150 ABIG = AB I G + DLOGIO(BIG) 

0036100 IF (ABIG.GT.EPSMAX) GOTO 55 

0036200 155 I F ( I COUNT. GE . 1 0 ) GOTO 55 

0036300 DO 54 1=1, NM1 

0036400 IF (.NOT. INTHd )) GOTO 54 

0036500 TEMP = VECT( I ) 

0036600 VECT ( I ) = VECT(I+1) 

0036700 VECT ( I + 1 ) = TEMP 

0036800 54 VECT(I+1) = VECT ( I +1 ) +MULT ( I ) *VECT ( I ) 

0036900 ICOUNT = I COUNT + 1 

0037000 GO TO 49 

0037100 55 IF (M.LE.2) GOTO 69 

0037200 MCOUNT(L) = ICOUNT 

0037300 MM2 = M-2 

0037400 DO 57 1=1, MM2 

0037500 Mil = M-l-l 

0037600 Mil = M-l+1 

0037700 DO 56 J=MI1,M 

0037800 56 VECT ( J ) =H ( J ,M1 I ) *VECT(M1 I +1 ) +VECT( J ) 

0037900 I NDEX = JNT(M1I ) 
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0038000 


TEMP = VECT (Ml I +1 ) 


0038100 


VECT(M1 1+1) = VECT( INDEX) 


0038200 

57 

VECT (INDEX) = TEMP 


0038300C 




0038400C 


NORMALIZE EIGENVECTOR. 


0038500C 




0038600 

69 

SUM = 0. 


0038700 


DO 58 l=l,M 


0038800 


RELVEC = VECT ( 1 ) 


0038900 


AMGVEC = CONJ 1 *VECT ( 1 ) 


0039000 

58 

SUM = SUM + RELVEC*RELVEC + AMGVEC*AMGVEC 

0039100 


SUM = DSQRT (SUM) 


0039200 


IF (SUM.EQ.O.) GO TO 60 


0039300 


DO 59 1=1, M 


0039400 

59 

VECT ( 1 ) = VECT ( 1 ) /SUM 


0039500 

60 

CONTINUE 


0039600 


DO 61 1 =1/ M 


0039700 

61 

A(I,L) = VECT ( 1 ) 


0039800 


CALL CPUTIM(KTIM) 


0039900 


KTSUM = KTSUM + (KTIM - JTIM) 


0040000 

63 

CONTINUE 


0040100 


NCAL = NCAL + NC 


0040200 

64 

IF(NSTART.EQ.l) GOTO 70 


0040300 


SHI FT ( 1 ) = 0. 


0040400 


NSTOP = NS TART - 1 


0040500 


N = NSTOP 


0040600 


GO TO 20 


0040700 

70 

DO 80 L=2/ M 


0040800 


DO 79 1 = 1, M 


0040900 

79 

JNT(I) = 0 


0041000 


RELAML = LAMBDA ( L ) 


0041100 


AM GAM L = CONJ 1 *LAMBDA( L) 


0041200 


LM1 = L - 1 


0041300 


R = 0 


0041400 


DO 71 1=1, LM1 


0041500 


RELAM 1 = LAMBDA ( 1 ) 


0041600 


AMGAMI = CONJ 1 *LAMBPA( 1 ) 


0041700 


IF ( DABS ( RELAML-RELAM 1 ) . GT . EPS ) GOTO 

71 

0041800 


IF (DABS ( AMGAML-AMGAM 1 ).GT. EPS) GOTO 

71 

0041900 


JNT { I ) = L 


0042000 


R = R + 1 


0042100 

71 

CONTINUE 


0042200 


IF (R.EQ.O) GOTO 80 


0042300 


DO 72 1=1, M 


0042400 

72 

VECT(I) = 0. 


0042500 


DO 75 1=1, LM1 


0042600 


IF (JNT ( 1 ).NE.L) GOTO 75 


0042700 


TEMP = 0. 


0042800 


DO 73 J=1,M 


0042900 

73 

TEMP = TEMP + DCONJG(A(J,L))*A(J, 1 ) 


0043000 


DO 74 J=1,M 


0043100 

74 

VECT(J) = VECT(J) + TEMP*A( J, 1 ) 


0043200 


IF (R.EQ.l) GOTO 76 


0043300 


R = R ^ 1 


0043400 

75 

CONTINUE 


0043500 

76 

SUM = 0. 
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0043600 


DO 77 1=1, M 

0043700 


A( 1, L) = A( 1 , L) - VECT ( 1 ) 

0043800 

77 

SUM = SUM + A(l,L)*DCONJG(A( 1,0) 

0043900 


IF (SUM.EQ.O.) GOTO 80 

0044000 


SUM = DSQRT ( SUM) 

0044100 


DO 78 1=1, M 

0044200 

78 

A( 1 , L ) = A ( 1 , L ) /SUM 

0044300 

80 

CONTINUE 

0044400 

92 

DO 95 J=l, M 

0044500 


DO 95 1=1, M 

0044600 


TEMP = AC 1 , J ) 

0044700 


AC 1 , J ) = AA( 1, J) 

0044800 

95 

AA ( 1 , J ) = TEMP 

0044900 


RETURN 

0045000 


ENTRY EVDATA ( 1 TS, KTS, NCO,MCO, RNORM) 

0045100 


DIMENSION MCO (l),NCO(l), RNORM (1) 

0045200 


DO 83 1=1, M 

0045300 

83 

NCO( 1 ) = NCOUNT ( 1 ) 

0045400 


ITS = ITSUM 

0045500 


IF ( . NOT. EVECT) RETURN 

0045600 


DO 84 1=1, M 

0045700 

84 

MCOC 1 ) = MCOUNT ( 1 ) 

0045800 


ANORM = 0. 

0045900 


DO 85 1=1, M 

0046000 


DO 85 J=1,M 

0046100 

85 

ANORM = ANORM + A( J, 1 ) *DCON JG ( A ( J, 1 ) ) 

0046200 


ANORM = DSQRT (ANORM ) 

0046300 


IF (ANORM. EQ.O.) ANORM=l. 

0046400 


KTS = KTSUM 

0046500 


DO 90 L=1,M 

0046600 


VNORM = 0. 

0046700 


DO 89 1=1, M 

0046800 


TEMP = 0. 

0046900 


DO 82 J = 1 , M 

0047000 

82 

TEMP = TEMP + A ( 1 , J ) *AA( J, L ) 

0047100 


TEMP = TEMP - LAMBDA ( L ) *AA ( 1 , L) 

0047200 

89 

VNORM = VNORM + ( CDABS (TEMP ) ) * * 2 

0047300 

90 

RNORM ( L ) = DSQRT ( VNORM ) /ANORM 

0047400 


RETURN 

0047500 


END 
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